home *** CD-ROM | disk | FTP | other *** search
/ Nebula 1 / Nebula One.iso / Graphics / Misc / NeXTcontour_1.7 / Source / contour_plot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-12  |  12.8 KB  |  649 lines
  1. /* contour_plot.f -- translated by f2c (version of 22 January 1991  19:25:10).
  2.    You must link the resulting object file with the libraries:
  3.     -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #import <appkit/color.h>
  7. #import "f2c.h"
  8. #import <dpsclient/wraps.h>
  9.  
  10.  
  11. /* Common Block Declarations */
  12.  
  13. struct {
  14.     real xt[1000], yt[1000], zt[1000];
  15.     integer ia[3000];
  16. } temp1_;
  17.  
  18. #define temp1_1 temp1_
  19.  
  20. /* Table of constant values */
  21.  
  22. static integer c__1 = 1;
  23.  
  24. /* Subroutine */ int contour_(jdim, kdim, nj, nk, x, y, f, ncont, acont,
  25.                   ppxunit, ppyunit, colorMap)
  26. integer *jdim, *kdim, *nj, *nk, *colorMap;
  27. real *x, *y, *f;
  28. integer *ncont;
  29. real *acont;
  30. float *ppxunit, *ppyunit;
  31. {
  32.     /* System generated locals */
  33.     integer x_dim1, x_offset, y_dim1, y_offset, f_dim1, f_offset;
  34.  
  35.     /* Local variables */
  36.     extern /* Subroutine */ int con2l_();
  37.  
  38.     /* Parameter adjustments */
  39.     --acont;
  40.     f_dim1 = *jdim;
  41.     f_offset = f_dim1 + 1;
  42.     f -= f_offset;
  43.     y_dim1 = *jdim;
  44.     y_offset = y_dim1 + 1;
  45.     y -= y_offset;
  46.     x_dim1 = *jdim;
  47.     x_offset = x_dim1 + 1;
  48.     x -= x_offset;
  49.  
  50.     /* Function Body */
  51.     con2l_(jdim, kdim, &c__1, nj, &c__1, nk, &x[x_offset], &y[y_offset], &f[
  52.         f_offset], ncont, &acont[1], ppxunit, ppyunit, colorMap);
  53.  
  54.     return 0;
  55. } /* contour_ */
  56.  
  57. /* Subroutine */ int cline2_(cont, np, x, y, 
  58.                  ppxunit, ppyunit, colorMap, icont, ncont)
  59. real *cont;
  60. integer *np, *colorMap;
  61. int icont, *ncont;
  62. real *x, *y;
  63. float *ppxunit, *ppyunit;
  64. {
  65.   int i;
  66.   float pattern0[] = {};    /* solid      */
  67.   float pattern1[] = {4.0, 4.0}; /* dash       */
  68.   NXColor color;
  69.   float r,g,b;
  70.   float h,s,br;
  71.   int one_third;
  72.  
  73.     /* Parameter adjustments */
  74.     --y;
  75.     --x;
  76.  
  77. /*  PSsetlinewidth(0.0);  set outside ! */
  78.  
  79. /*  PSsetgray(0);  set outside ! */
  80.  
  81.   switch(*colorMap) {
  82.   case 0:
  83.     PSsetdash(pattern0, 0, 0.0);
  84.     break;
  85.   case 1:
  86.     if(*cont >= 0.0)PSsetdash(pattern0, 0, 0.0);
  87.     if(*cont < 0.0) PSsetdash(pattern1, 2, 0.0);
  88.     break;
  89.   case 2:
  90.     PSsetdash(pattern0, 0, 0.0);
  91.     one_third = *ncont/3;
  92.     r = 0;
  93.     g = 0;
  94.     b = 0;
  95.     h = 0;
  96.     s = 0;
  97.     br = 0;
  98. /*    if(icont <= one_third)
  99.       {
  100.     b = 1.-(float)icont/one_third;
  101.     g = (float)icont/one_third;
  102.     r = b*g;
  103.       }
  104.     if(icont > one_third && icont <= 2*one_third)
  105.       {
  106.     b = ((float)icont - (float)one_third)/one_third;
  107.     g = 1. - ((float)icont - (float)one_third)/one_third;
  108.     r = g*b;
  109.       }
  110.     if(icont > 2*one_third && icont <= *ncont)
  111.       {
  112.     r = 1.-((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
  113.     b = ((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
  114.     g = r*b;
  115.       }
  116.  
  117.     color = NXConvertRGBToColor(r,g,b);
  118. */
  119.     h = 0.6666*icont/(float)*ncont;
  120.     s = 1.0;
  121.     br = 1.0;
  122.     color = NXConvertHSBToColor(h,s,br);
  123.        
  124.     NXSetColor(color);
  125.     break;
  126.  
  127.   case 3:
  128.     if(*cont >= 0.0)PSsetdash(pattern0, 0, 0.0);
  129.     if(*cont < 0.0) PSsetdash(pattern1, 2, 0.0);
  130.     one_third = *ncont/3;
  131.     r = 0;
  132.     g = 0;
  133.     b = 0;
  134.     h = 0;
  135.     s = 0;
  136.     br = 0;
  137. /*    if(icont <= one_third)
  138.       {
  139.     b = 1.-(float)icont/one_third;
  140.     g = (float)icont/one_third;
  141.     r = b*g;
  142.       }
  143.     if(icont > one_third && icont <= 2*one_third)
  144.       {
  145.     b = ((float)icont - (float)one_third)/one_third;
  146.     g = 1. - ((float)icont - (float)one_third)/one_third;
  147.     r = g*b;
  148.       }
  149.     if(icont > 2*one_third && icont <= *ncont)
  150.       {
  151.     r = 1.-((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
  152.     b = ((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
  153.     g = r*b;
  154.       }
  155.  
  156.     color = NXConvertRGBToColor(r,g,b);
  157. */
  158.     h = 0.6666*icont/(float)*ncont;
  159.     s = 1.0;
  160.     br = 1.0;
  161.     color = NXConvertHSBToColor(h,s,br);
  162.        
  163.     NXSetColor(color);
  164.     break;
  165.   }
  166.  
  167.   PSnewpath();
  168.   PSmoveto(x[1]*(*ppxunit), y[1]*(*ppyunit));
  169.   for (i=2; i <= *np; i++){
  170.     PSlineto(x[i]*(*ppxunit), y[i]*(*ppyunit));
  171.   }
  172.   PSstroke();
  173.  
  174.  
  175.     /* Function Body */
  176.     return 0;
  177. } /* cline2_ */
  178.  
  179.  
  180.  
  181. /* ----------------------------------------------------------------------- */
  182. /* Subroutine */ int con2l_(idim, jdim, is, ie, js, je, x, y, f, ncont, acont,
  183.                 ppxunit, ppyunit, colorMap)
  184.  
  185. integer *idim, *jdim, *is, *ie, *js, *je, *colorMap;
  186. real *x, *y, *f;
  187. integer *ncont;
  188. real *acont;
  189. float *ppxunit, *ppyunit;
  190. {
  191.     /* System generated locals */
  192.     integer x_dim1, x_offset, y_dim1, y_offset, f_dim1, f_offset, i__1, i__2, 
  193.         i__3, i__4;
  194.  
  195.     /* Local variables */
  196.     static integer ibeg, jbeg, idir;
  197.     static real cont;
  198.     static integer i, j, iaend, icont;
  199.     extern /* Subroutine */ int cline2_();
  200.     static integer ij, np;
  201.     static real wx, wy;
  202.     static integer ignext, lnstrt, iae, iia, iee, jee, imb, ima;
  203.     static real fij;
  204.     static integer iss, jss;
  205.     static real xyf, wxx, wyy;
  206.  
  207.  
  208. /*  Calculate contour lines for the function F in the region IS to IE, JS 
  209. to*/
  210. /*  JE.  X,Y coordinates corresponding to the grid points are in arrays X 
  211. and*/
  212. /*   Y. */
  213.  
  214. /*  Common block TEMP1 contains scratch arrays XT, YT, ZT, and IA, and may
  215.  be*/
  216. /*  used elsewhere.  This subroutine is supposed to be able to recover fro
  217. m*/
  218. /*   filling either array. */
  219.  
  220.  
  221. /*   One little check.  If IS=IE or JS=JE, return with no contour lines. 
  222. */
  223.  
  224.     /* Parameter adjustments */
  225.     --acont;
  226.     f_dim1 = *idim;
  227.     f_offset = f_dim1 + 1;
  228.     f -= f_offset;
  229.     y_dim1 = *idim;
  230.     y_offset = y_dim1 + 1;
  231.     y -= y_offset;
  232.     x_dim1 = *idim;
  233.     x_offset = x_dim1 + 1;
  234.     x -= x_offset;
  235.  
  236.     /* Function Body */
  237.     if (*is == *ie || *js == *je) {
  238.     return 0;
  239.     }
  240.  
  241. /*   Zero the marker array. */
  242.  
  243.     for (iia = 1; iia <= 3000; ++iia) {
  244.     temp1_1.ia[iia - 1] = 0;
  245. /* L100: */
  246.     }
  247.  
  248. /*   LNSTRT=1 means we're starting a new line. */
  249.  
  250.     lnstrt = 1;
  251.     ignext = 0;
  252.  
  253. /*   Loop through each contour level. */
  254.  
  255.     i__1 = *ncont;
  256.     for (icont = 1; icont <= i__1; ++icont) {
  257.     cont = acont[icont];
  258.  
  259. /*  The IS-IE,JS-JE region may have to be subdivided because of IA spa
  260. ce.  This*/
  261. /*  may have a detrimental effect on labelling of the contour lines.  
  262. Move*/
  263. /*   IS,IE,JS,JE to new variables. */
  264.  
  265.     iss = *is;
  266.     iee = *ie;
  267.     jss = *js;
  268.     jee = *je;
  269.  
  270. /*  Flag points in IA where the the function increases through the con
  271. tour*/
  272. /*  level, not including the boundaries.  This is so we have a list of
  273.  at least*/
  274. /*   one point on each contour line that doesn't intersect a boundary.
  275.  */
  276.  
  277. L110:
  278.     iae = 0;
  279.     i__2 = jee - 1;
  280.     for (j = jss + 1; j <= i__2; ++j) {
  281.         imb = 0;
  282.         iaend = iae;
  283.         i__3 = iee;
  284.         for (i = iss; i <= i__3; ++i) {
  285.         if (f[i + j * f_dim1] <= cont) {
  286.             imb = 1;
  287.         } else if (imb == 1) {
  288.             ++iae;
  289.             temp1_1.ia[iae - 1] = i * 1000 + j;
  290.             imb = 0;
  291.  
  292. /*  Check if the IA array is full.  If so, the subdividing
  293.  algorithm goes like*/
  294. /*  this:  if we've marked at least one J row, drop back t
  295. o the last completed*/
  296. /*  J and call that the region.  If we haven't even finish
  297. ed one J row, our*/
  298. /*   region just extends to this I location. */
  299.  
  300.             if (iae == 3000) {
  301.             if (j > jss + 1) {
  302.                 iae = iaend;
  303.                 jee = j;
  304.             } else {
  305. /* Computing MIN */
  306.                 i__4 = j + 1;
  307.                 jee = min(i__4,jee);
  308.                 iee = i;
  309.             }
  310.             goto L210;
  311.             }
  312.  
  313.         }
  314. /* L200: */
  315.         }
  316.     }
  317.  
  318. /*   Search along the boundaries for contour line starts.  IMA tells w
  319. hich */
  320. /*   boundary of the region we're on. */
  321.  
  322. L210:
  323.     ima = 1;
  324.     imb = 0;
  325.     ibeg = iss - 1;
  326.     jbeg = jss;
  327.  
  328. L1:
  329.     ++ibeg;
  330.     if (ibeg == iee) {
  331.         ima = 2;
  332.     }
  333.     goto L5;
  334. L2:
  335.     ++jbeg;
  336.     if (jbeg == jee) {
  337.         ima = 3;
  338.     }
  339.     goto L5;
  340. L3:
  341.     --ibeg;
  342.     if (ibeg == iss) {
  343.         ima = 4;
  344.     }
  345.     goto L5;
  346. L4:
  347.     --jbeg;
  348.     if (jbeg == jss) {
  349.         ima = 5;
  350.     }
  351. L5:
  352.     if (f[ibeg + jbeg * f_dim1] > cont) {
  353.         goto L7;
  354.     }
  355.     imb = 1;
  356. L6:
  357.     switch ((int)ima) {
  358.         case 1:  goto L1;
  359.         case 2:  goto L2;
  360.         case 3:  goto L3;
  361.         case 4:  goto L4;
  362.         case 5:  goto L91;
  363.     }
  364. L7:
  365.     if (imb != 1) {
  366.         goto L6;
  367.     }
  368.  
  369. /*   Got a start point. */
  370.  
  371.     imb = 0;
  372.     i = ibeg;
  373.     j = jbeg;
  374.     fij = f[ibeg + jbeg * f_dim1];
  375.  
  376. /*   Round the corner if necessary. */
  377.  
  378.     switch ((int)ima) {
  379.         case 1:  goto L21;
  380.         case 2:  goto L11;
  381.         case 3:  goto L12;
  382.         case 4:  goto L13;
  383.         case 5:  goto L51;
  384.     }
  385. L11:
  386.     if (j != jss) {
  387.         goto L31;
  388.     }
  389.     goto L21;
  390. L12:
  391.     if (i != iee) {
  392.         goto L41;
  393.     }
  394.     goto L31;
  395. L13:
  396.     if (j != jee) {
  397.         goto L51;
  398.     }
  399.     goto L41;
  400.  
  401. /*   This is how we start in the middle of the region, using IA. */
  402.  
  403. L10:
  404.     i = temp1_1.ia[iia - 1] / 1000;
  405.     j = temp1_1.ia[iia - 1] - i * 1000;
  406.     fij = f[i + j * f_dim1];
  407.     temp1_1.ia[iia - 1] = 0;
  408.     goto L21;
  409. /*                                                                    
  410.    4 */
  411. /*  Look different directions to see which way the contour line went: 
  412. 1-|-3*/
  413. /*                                                                    
  414.    2 */
  415. L20:
  416.     ++j;
  417. L21:
  418.     --i;
  419.     if (i < iss) {
  420.         goto L90;
  421.     }
  422.     idir = 1;
  423.     if (f[i + j * f_dim1] <= cont) {
  424.         goto L52;
  425.     }
  426.     fij = f[i + j * f_dim1];
  427.     goto L31;
  428. L30:
  429.     --i;
  430. L31:
  431.     --j;
  432.     if (j < jss) {
  433.         goto L90;
  434.     }
  435.     idir = 2;
  436.     if (f[i + j * f_dim1] <= cont) {
  437.         goto L60;
  438.     }
  439.     fij = f[i + j * f_dim1];
  440.     goto L41;
  441. L40:
  442.     --j;
  443. L41:
  444.     ++i;
  445.     if (i > iee) {
  446.         goto L90;
  447.     }
  448.     idir = 3;
  449.     if (f[i + j * f_dim1] <= cont) {
  450.         goto L60;
  451.     }
  452.     fij = f[i + j * f_dim1];
  453.     goto L51;
  454. L50:
  455.     ++i;
  456. L51:
  457.     ++j;
  458.     idir = 4;
  459.     if (j > jee) {
  460.         goto L90;
  461.     }
  462.     if (f[i + j * f_dim1] <= cont) {
  463.         goto L60;
  464.     }
  465.     fij = f[i + j * f_dim1];
  466.     goto L21;
  467.  
  468. /*   Wipe this point out of IA if it's in the list. */
  469.  
  470. L52:
  471.     if (iae == 0) {
  472.         goto L60;
  473.     }
  474.     ij = i * 1000 + j + 1000;
  475.     i__3 = iae;
  476.     for (iia = 1; iia <= i__3; ++iia) {
  477.         if (temp1_1.ia[iia - 1] == ij) {
  478.         temp1_1.ia[iia - 1] = 0;
  479.         goto L60;
  480.         }
  481. /* L300: */
  482.     }
  483.  
  484. /*   Do interpolation for X,Y coordinates. */
  485.  
  486. L60:
  487.     xyf = (cont - f[i + j * f_dim1]) / (fij - f[i + j * f_dim1]);
  488.  
  489. /*  This tests for a contour point coinciding with a grid point.  In t
  490. his case*/
  491. /*  the contour routine comes up with the same physical coordinate twi
  492. ce.  If*/
  493. /*  If we don't trap it, it can (in some cases significantly) increase
  494.  the*/
  495. /*  number of points in a contour line.  Also, if this happens on the 
  496. first*/
  497. /*  point in a line, the second point could be misinterpreted as the e
  498. nd of a*/
  499. /*   (circling) contour line. */
  500.  
  501.     if (xyf == (float)0.) {
  502.         ++ignext;
  503.     }
  504.     switch ((int)idir) {
  505.         case 1:  goto L61;
  506.         case 2:  goto L62;
  507.         case 3:  goto L63;
  508.         case 4:  goto L64;
  509.     }
  510. L61:
  511.     wxx = x[i + j * x_dim1] + xyf * (x[i + 1 + j * x_dim1] - x[i + j * 
  512.         x_dim1]);
  513.     wyy = y[i + j * y_dim1] + xyf * (y[i + 1 + j * y_dim1] - y[i + j * 
  514.         y_dim1]);
  515.     goto L65;
  516. L62:
  517.     wxx = x[i + j * x_dim1] + xyf * (x[i + (j + 1) * x_dim1] - x[i + j * 
  518.         x_dim1]);
  519.     wyy = y[i + j * y_dim1] + xyf * (y[i + (j + 1) * y_dim1] - y[i + j * 
  520.         y_dim1]);
  521.     goto L65;
  522. L63:
  523.     wxx = x[i + j * x_dim1] + xyf * (x[i - 1 + j * x_dim1] - x[i + j * 
  524.         x_dim1]);
  525.     wyy = y[i + j * y_dim1] + xyf * (y[i - 1 + j * y_dim1] - y[i + j * 
  526.         y_dim1]);
  527.     goto L65;
  528. L64:
  529.     wxx = x[i + j * x_dim1] + xyf * (x[i + (j - 1) * x_dim1] - x[i + j * 
  530.         x_dim1]);
  531.     wyy = y[i + j * y_dim1] + xyf * (y[i + (j - 1) * y_dim1] - y[i + j * 
  532.         y_dim1]);
  533.  
  534. /*   Figure out what to do with this point. */
  535.  
  536. L65:
  537.     if (lnstrt != 1) {
  538.         goto L66;
  539.     }
  540.  
  541. /*   This is the first point in a contour line. */
  542.  
  543.     np = 1;
  544.     temp1_1.xt[np - 1] = wxx;
  545.     temp1_1.yt[np - 1] = wyy;
  546.     wx = wxx;
  547.     wy = wyy;
  548.     lnstrt = 0;
  549.     goto L67;
  550.  
  551. /*   Add a point to this line.  Check for duplicate point first. */
  552.  
  553. L66:
  554.     if (ignext == 2) {
  555.         if (wxx == temp1_1.xt[np - 1] && wyy == temp1_1.yt[np - 1]) {
  556.         ignext = 0;
  557.         goto L67;
  558.         } else {
  559.         ignext = 1;
  560.         }
  561.     }
  562.  
  563.     ++np;
  564.     temp1_1.xt[np - 1] = wxx;
  565.     temp1_1.yt[np - 1] = wyy;
  566.  
  567. /*   See if the temporary array xT is full. */
  568.  
  569.     if (np == 1000) {
  570.         cline2_(&cont, &np, temp1_1.xt, temp1_1.yt, 
  571.             ppxunit, ppyunit, colorMap, icont, ncont);
  572.         temp1_1.xt[0] = temp1_1.xt[np - 1];
  573.         temp1_1.yt[0] = temp1_1.yt[np - 1];
  574.         np = 1;
  575.     }
  576.  
  577. /*   Check to see if we're back to the initial point. */
  578.  
  579.     if (wxx != wx) {
  580.         goto L67;
  581.     }
  582.     if (wyy == wy) {
  583.         goto L90;
  584.     }
  585.  
  586. /*   Nope.  Search for the next point on this line. */
  587.  
  588. L67:
  589.     switch ((int)idir) {
  590.         case 1:  goto L50;
  591.         case 2:  goto L20;
  592.         case 3:  goto L30;
  593.         case 4:  goto L40;
  594.     }
  595.  
  596. /*  This is the end of a contour line.  After this we'll start a new l
  597. ine.*/
  598.  
  599. L90:
  600.     lnstrt = 1;
  601.     ignext = 0;
  602.     cline2_(&cont, &np, temp1_1.xt, temp1_1.yt, 
  603.         ppxunit, ppyunit, colorMap, icont, ncont);
  604.  
  605. /*  If we're not done looking along the boundaries, go look there some
  606.  more.*/
  607.  
  608.     if (ima != 5) {
  609.         goto L6;
  610.     }
  611.  
  612. /*   Otherwise, get the next start out of IA. */
  613.  
  614. L91:
  615.     if (iae != 0) {
  616.         i__3 = iae;
  617.         for (iia = 1; iia <= i__3; ++iia) {
  618.         if (temp1_1.ia[iia - 1] != 0) {
  619.             goto L10;
  620.         }
  621. /* L400: */
  622.         }
  623.     }
  624.  
  625. /*  And if there are no more of these, we're done with this region.  I
  626. f we've*/
  627. /*   subdivided, update the region pointers and go back for more. */
  628.  
  629.     if (iee == *ie) {
  630.         if (jee != *je) {
  631.         jss = jee;
  632.         jee = *je;
  633.         goto L110;
  634.         }
  635.     } else {
  636.         iss = iee;
  637.         iee = *ie;
  638.         goto L110;
  639.     }
  640.  
  641. /*   Loop back for the next contour level. */
  642.  
  643. /* L1000: */
  644.     }
  645.  
  646.     return 0;
  647. } /* con2l_ */
  648.  
  649.